#############################################
#############################################
###############   Problem 5   ############### 


#############################################
############## Part (a) #####################

## Scan in data ##
y1 = c(2.34, 2.46, 2.83, 2.04, 2.69);
y2 = c(2.64, 3.00, 3.19, 3.83);
y3 = c(2.61, 2.07, 2.80, 2.58, 2.98, 2.30);
y4 = c(1.32, 1.62, 1.92, 0.88, 1.50, 1.30);
y5 = c(0.41, 0.83, 0.53, 0.32, 1.62);

## Create list of data to make calculations easier ##
data = list(y1, y2, y3, y4, y5);

## To obtain degress of freedom ##
t = 5;  ## Number of treatments
N1 = length(y1);
N2 = length(y2);
N3 = length(y3);
N4 = length(y4);
N5 = length(y5);
N = c(N1, N2, N3, N4, N5);

n = sum(N); ## total number of observations

## Calculate grand mean for later ##
g.mean = (sum(y1) + sum(y2) + sum(y3) + sum(y4) + sum(y5)) / n;

SSR = SSE = 0;

## Sums of Squares ##
for(i in 1:t){

	SSR = SSR + N[i] * (mean(data[[i]]) - g.mean)^2;
	SSE = SSE + sum((data[[i]] - mean(data[[i]]))^2);
}

SST = SSR + SSE;

## Calculate mean squares ##
MSR = SSR / (t - 1);
MSE = SSE / (n - t);

## Calculate F stat ##
F = MSR / MSE;

## Calculate the critical value ##
f = qf(1 - 0.05, t - 1, n - t);


#############################################
############## Part (b) #####################

lambda = c(1/3, 1/3, 1/3, -1/2, -1/2);

tmp = 0;
tmp2 = 0;

for(i in 1:t){

	tmp = tmp + lambda[i] * mean(data[[i]]);
	tmp2 = tmp2 + lambda[i]^2 / N[i];

}

F = (tmp^2 / tmp2) / MSE;
f = qf(1 - 0.05, 1, n - t);


#############################################
############## Part (c) #####################

## Calculate the sums of squares for contrasts ##
## First contrast ##
lambda = c(1/3, 1/3, 1/3, -1/2, -1/2);

tmp = 0;
tmp2 = 0;

for(i in 1:t){

	tmp = tmp + lambda[i] * mean(data[[i]]);
	tmp2 = tmp2 + lambda[i]^2 / N[i];

}

SS = (tmp^2 / tmp2);

lambda1 = c(0, 1, -1, 0, 0);

tmp1 = 0;
tmp12 = 0;

for(i in 1:t){

	tmp1 = tmp1 + lambda1[i] * mean(data[[i]]);
	tmp12 = tmp12 + lambda1[i]^2 / N[i];

}

SS1 = (tmp1^2 / tmp12);

## Next contrast ##
lambda2 = c(0, 0, 0, 1, -1);

tmp2 = 0;
tmp22 = 0;

for(i in 1:t){

	tmp2 = tmp2 + lambda2[i] * mean(data[[i]]);
	tmp22 = tmp22 + lambda2[i]^2 / N[i];

}

SS2 = (tmp2^2 / tmp22);

## Next contrast ##
lambda3 = c(1, -1/2, -1/2, 0, 0);

tmp3 = 0;
tmp32 = 0;

for(i in 1:t){

	tmp3 = tmp3 + lambda3[i] * mean(data[[i]]);
	tmp32 = tmp32 + lambda3[i]^2 / N[i];

}

SS3 = (tmp3^2 / tmp32);


#############################################
############## Part (d) #####################

lambda = c(1,-1);
N = c(6, 5);

LB = (mean(data[[4]]) - mean(data[[5]])) - qt(.975, 21) * sqrt(MSE * sum(lambda^2 / N));
UB = (mean(data[[4]]) - mean(data[[5]])) + qt(.975, 21) * sqrt(MSE * sum(lambda^2 / N));